# save outputs of array with variations 
# in penalty lambda, alternative Bn, dimension p, 

###################################################################
# MC 
###################################################################

rm(list = ls(all = T));
start_T <- Sys.time() ; 

library(foreach); library(doParallel)

###- controls -###
G = c(5*10^4) # of grid points for gamma
p0 = c(2, 5) # dimension of IV
cltype = "FORK"

###################

n = 1000  # sample size
B = 5000   # of iterations to compute power, H0:5000, H1:1000

bt0 = 1;  
Bn0 = c(0,.2)  # true parameter values

grid_ld = seq(0,1,length.out = 11) # penalty 

colP <- vector() ; 

for (j_p in 1:length(p0)) {   
  p = p0[j_p]

  all_T = array(NA,dim = c(length(Bn0),length(grid_ld),B)) 
  filename = paste("Power",p,(G/1000),".RData",sep = "_");

  for (H1 in 1:length(Bn0)) { 
  
    Bn = Bn0[H1]
    bt_n = bt0 + Bn
    
    numCores <- detectCores() - 3
    cl <- makeCluster(numCores,type = cltype) 
    clusterExport(cl, ls()) # export everything
    
    colT <- foreach(b = 1:B, .combine = 'cbind') %dopar% {
    
      grid_gm = matrix(runif(p*p*G,min=-5,max=5),p,G*p)
      grid_gm = cbind(grid_gm,-grid_gm)
      gm_1 = apply(abs(grid_gm), 2, sum)
      
      xx = matrix(rnorm(2*n),n,2); 
      ep13 = matrix(runif(2*n,min=-1,max = 1),n,2)
    
      x = xx[,1]/2 + ep13[,1];
      e = xx[,1]/2 + xx[,2]/2
      Y = x*bt_n + e
      U = Y - x*bt0;
    
      W = cbind(ep13[,1] - ep13[,2], matrix(runif(n*(p-1),min = -1,max = 1), n ,p-1))
      eWgm = exp(W %*% grid_gm) ;# dim(eWgm)
    
      Mn = apply(U*eWgm,2,mean)
      sn = sqrt(apply((U*eWgm)^2,2,mean))
      Qn = sqrt(n)*abs(Mn)/sn
    
      Tn = vector()
      for (j in 1:length(grid_ld)){
        Tn = c(Tn,max(Qn - grid_ld[j]*gm_1))
        }
    
      Tn
      
    }
    
    stopCluster(cl)
    all_T[H1,,] <- colT
    
    if (Bn == 0) {
      cv = apply(colT, 1, quantile,probs=0.95)
      saveRDS(cv,file= paste("cv","_W",p,".rds",sep = "") )
      
      }else{
    
      cv <- readRDS(paste("cv","_W",p,".rds",sep = ""))
      pwr = apply((colT > cv), 1, mean)
      colP <- cbind(colP,pwr)
      
      plot(grid_ld,pwr)
      }
    
    end_T <- Sys.time()
    elapsed_T <- end_T - start_T 
    
    save.image(paste("./",filename,sep=""))
    
  }
  
  difftime(end_T , start_T, units = "hour")
  
}

difftime(end_T , start_T, units = "hour")


######################################################################
# PLOT 
######################################################################
filename <- "Power_5_50_.RData"
load(paste("./",filename,sep=""))

Bn = 0.1;
pwrB = 1:2

hght = 4.5;
wdth = 6
test_col <- c(1:length(p0))
test_lin <- 1:length(p0) # c(3,2,1,4:n_t)
test_pch <- 15:19 #point symbols
cexleg <- 1.25

b_main = paste("Power function for each p: "," n =",n,sep = " ")

xlab = "penaly (10 * lambda) " 
ylab = "rejection rate"

pdf(paste("Figure_1",".pdf",sep="_"), height = hght, width = wdth)
plot(ts(colP[,pwrB],start = 0, end = 10, frequency = 1),
     type = "b", plot.type = "single", 
     pch = test_pch,
     col=test_col,lty=test_lin,lwd=2,cex = 1,
     main = b_main,cex.main=1, 
     xlab = xlab, ylab =ylab,cex.lab = 1,
     ylim = c(0.25,0.6)
    )
#abline(h=0,lty=3)
legend("topright",
       legend=paste("p=",as.vector(p0),sep = ""),
       col=test_col, pch = test_pch,cex=cexleg)
dev.off()
